paquetes <- c("GGally", "HistData", "RColorBrewer", "caret", "CDR", "classInt", "corrplot","data.table", "DiagrammeR", "devtools", "DT", "eurostat", "explore", "factoextra", "gapminder", "ggalluvial", "gganimate", "ggmosaic", "ggplot2", "ggpubr", "ggrepel", "ggridges", "ggstatsplot", "gplots", "gstat", "ggthemes", "igraph", "knitr", "leaflet", "leaflet.extras", "leaflet.providers", "lmtest", "mapSpain", "openxlsx", "pacman", "pROC", "plotly", "psych", "quarto", "raster", "readxl", "remotes", "rlang", "rvest", "reshape2", "scales", "sf", "sjPlot","skimr", "summarytools", "survival", "tidyverse", "titanic", "tm", "visdat", "viridis", "wordcloud", "quantmod","fBasics","curl","dygraphs","tidyquant", "timetk","highcharter", "wooldridge")
instalarpaquetes <- function(x){ # Creamos nuestra primera función :)
for (i in x){ # Para cada ítem en el objeto x
if (!require(i, character.only = TRUE)){
install.packages(i, dependencies = TRUE)
library(i, character.only = TRUE)
}
}
}
instalarpaquetes(paquetes)Introducción a R para la Economía y el Análisis de Políticas Públicas
Desde La Mancha con amor
0. Manuales interesantes
En este apartado os dejo información sobre algunos manuales y recursos que han sido utilizados para el desarrollo de este documento, así como otros que me han acompañado durante los últimos años.
Introducción a la programación en R
Libros de econometría y estadística aplicada
- Econometría aplicada utilizando R
- Algo más avanzado: R for HTA (no es libro como tal, pero tiene recursos)
- Baio et al. (2017). Bayesian Cost-Effectiveness Analysis with the R package BCEA.
- Estadística aplicada avanzada con R (buscar libro en lugares prohibidos)
- Cook y Lawless (2018). Multistate Models for the Analysis of Life History Data - Este está muy chulo para entender modelos multiestado que son bastante interesantes para entender transiciones entre estados (ej.: empleo, desempleo, etc.).
Gráficos y visualización
- ggplot2: elegant graphics for Data Analysis - Utilizado para preparar este taller
Epidemiología para salud pública
- The Epidemiologist R Handbook: R for applied epidemiology and public health (Incluye Data Management, que suele dar problemas) - Utilizado para preparar este taller
Microeconometría y análisis económico
- Adams, C. (2021). Learning Microeconometrics with R - Utilizado para preparar este taller
- López Cano (2019). Análisis de datos con R aplicados a la economía, la empresa y la industria
Bioestadística y ciencia de datos
- Martínez González (2020). Bioestadística amigable. (Usa Stata y SPSS, pero sirve para conceptos)
- Fernández y Montero (2024). Fundamentos de Ciencia de Datos con R - Utilizado para preparar este taller
- Zimmer, Powell y Velásquez (2024). Exploring Complex Survey Data Analysis using R
- Chan (2016). Biostatistics for Epidemiology and Public Health Using R
Otros
- R for Data Science Este es un clásico joven que han hecho los creadores del tidyverse y también tienen cositas en ggplot2 - gente fiera, vaya.
- Advanced R
- Hands-On Programming with R
1. Recordatorio de instalación de R, RStudio y Quarto
- Instalación de R y RStudio: https://posit.co/download/rstudio-desktop/
- Instalar Rtools y devtools: Para descargar paquetes de github
- Devtools: install.packages(“devtools”)
- Rtools: https://cloud.r-project.org/bin/windows/Rtools/rtools44/rtools.html /
- Instalar Quarto: https://quarto.org/docs/get-started/
Paquetes que (quizá) utilizaremos
Antes de avanzar, recuerda es relevante fijar el directorio de trabajo. Si estás en un proyecto de R o de Quarto el directorio se ubica en ese mismo proyecto. Si simplemente trabajas con carpetas en Windows (no tiene por qué estar mal hacerlo así), recuerda fijarlo en la carpeta correspondiente.Se puede predefinir en las Global Options de RStudio.
getwd() # Muestra el directorio de trabajo[1] "C:/Users/Roberto.MLacoba/OneDrive - Universidad de Castilla-La Mancha/UCLM/Clases universitarias/Lo básico en R/Intro R Eco y APP"
#setwd() # Fija el directorio de trabajo en la ruta asignada, por ejemplo: C:/Users/Usuario/Desktop2. Operadores elementales e iniciación a la programación
2.1 Operadores elementales
2.1.1 Operadores aritméticos
- Suma:
+ - Resta:
- - Multiplicación:
* - División:
/ - Exponenciación:
^ - Módulo:
%%(resto de la división) - División entera:
%/%(cociente de la división) - Raíz cuadrada:
sqrt() - Logaritmo natural:
log() - Logaritmo en base 10:
log10() - Exponencial:
exp()e^x - Valor absoluto:
abs() - Factorial:
factorial() - Producto de matrices:
%*%
# Ejemplos de operadores aritméticos
# Ejemplos de operadores aritméticos
x <- 5
y <- 3
z <- 2
x+y # Suma[1] 8
x-y # Resta[1] 2
sqrt(y)[1] 1.732051
exp(1)[1] 2.718282
log(10)[1] 2.302585
log10(100)[1] 2
[1] 1
a<- as.matrix(c(1,2,3)) #3x1
b<- as.matrix(c(4,5,6)) #3x1
bt<-t(as.matrix(c(4,5,6))) #1x3)
a%*%b # No funcionaError in a %*% b: argumentos no compatibles
a%*%bt # 3x1 * 1x3 = 3x3 [,1] [,2] [,3]
[1,] 4 5 6
[2,] 8 10 12
[3,] 12 15 18
bt%*%a # 1x3 * 3x1 = 1x1 [,1]
[1,] 32
2.1.2 Operadores de comparación
- Mayor que:
> - Menor que:
< - Mayor o igual que:
>= - Menor o igual que:
<= - Igual que:
== - Distinto de:
!=
2.1.3 Operadores lógicos
2.1.4 Operadores de asignación
- Asignación simple:
<-o=o-> - El pipe
%>%(de magrittr) o también|>(de base R): se puede usar para encadenar funciones y se interpreta como un “entonces”. Traducción literal:tubería
Vemos unos ejemplos con 2.1.2, 2.1.3 y 2.1.4 utilizando R base.
# Usamos la clásica base de datos iris
a<- mean(iris$Sepal.Length[iris$Sepal.Length>=5]) # Media de la longitud del sépalo si es mayor o igual que 5
a[1] 6.041406
b=mean(iris$Sepal.Length[iris$Sepal.Length>=5 & iris$Species=="setosa"]) # Media de la longitud del sépalo si es mayor o igual que 5 y la especie es setosa
b[1] 5.23
mean(iris$Sepal.Length[iris$Sepal.Length>=5 | iris$Petal.Length>=1.5]) ->c # Media de la longitud del sépalo si es mayor o igual que 5 o la longitud del pétalo es mayor o igual que 1.5
c[1] 5.960584
Como se observa, R ofrece una versatilidad razonable ya que permite subdividir de forma sencilla cualquier muestra/información. Imagina cómo harías esto en Excel… Es más, vamos a hacerlo y comparamos (y de paso aprendemos a exportar documentos en formato xlsx).
# Exportamos a Excel
library(openxlsx)
write.xlsx(iris, "data/iris.xlsx")2.1.5 Operadores de concatenación
[1] "Hola Roger: tot bé?"
2.2 Iniciación a la programación
Programar es un arte, un arte divertido, pero nadie nace aprendido y el ensayo y el error es clave (esto también aplica al manejo de R). En este apartado se recogen algunos elementos fundamentales para que sirva de introducción: tómalo como si fuese un juego. Puedes intentar reproducir el código a mano (así practicas escribiendo a mano) o cópialo, úsalo y altéralo para ver qué pasa. Lo que sí es un aspecto relevante de la programación (y no solo de la programación, también en otras disciplinas) es la posibilidad de imaginar. Cuando queramos hacer algún tipo de alteración de una base de datos debemos pensar qué es lo que realmente queremos hacer: ¿queremos obtener algo que responde de condiciones?, ¿queremos reproducir y repetir una serie de valores siguiendo un determinado patrón?, ¿queremos que ocurra algo si se cumple una condición? Pensemos, en este sentido, como con el lenguaje natural. Decía un filósofo, no recuerdo el nombre, que el mundo de las personas es tan amplio como las palabras que conocen. En estas disciplinas pensemos de forma similar: podremos hacer tanto como podamos ser capaces de imágenes; seguir formándonos y aprendiendo es fundamental. El diablo sabe más por viejo que por diablo.
2.2.1 Estructuras de control
- Condicionales:
if,else,else if - Bucles:
for,while,repeat - Interrupción de bucles:
break,next - Funciones:
function() - Estructuras de control:
switch(),tryCatch(). Esto sirve para evaluar una opción(switch)o para capturar errores(tryCatch).
Profundizaremos en los condicionales, en los bucles y en su interrupción.
2.2.2 Ejemplos de condicionales, bucles y de interrupción de bucles
Partimos de unos ejemplos sencillos sobre la propia consola/script de R con el fin de familiarizarnos con R y su entorno:
# Ejemplo de condicional
x <- 5
if (x > 5){
print("x es mayor que 5")
} else {
print("x es menor o igual que 5")
}
# Ejemplo de bucle
for (i in 1:5){
print(i)
}
# Ejemplo de bucle creando un vector:
vector <- c() # Definimos un vector vacío
for (i in rep(1:10,times=2, each=5)){
vector <- c(vector, i) # Lo que hace es iterar y añadir el valor de i al vector
}
print(vector)
vector <- c()
for (i in seq(1,10,2)){
vector <- c(vector, i)
}
print(vector)
# Ahora (muchísimo) más eficiente:
vector <- rep(1:10, times=2, each=5)
# Nota: pensar siempre cuál es la mejor forma de llegar al resultado
# Ejemplo de bucle con interrupción
for (i in 1:10){
if (i == 5){
break
}
print(i)
}
# Ejemplo con repeat:
i <- 1 # Número de partida
repeat{ # Repite
print(i) # Imprime i
i <- i + 1 # Suma 1 a i
if (i > 100){ # Si i es mayor que 100
break # Para
}
}
# Ejemplo con while:
i <- 1
while (i <= 10){ # Mientras i sea menor o igual que 10
print(i) # Imprime i
i <- i + 1 # Suma 1 a i
}
# Ejemplo combinando bucles y condicionales
for (i in 1:10){
if (i %% 2 == 0){ # Si el resto de dividir i entre 2 es 0
print(paste(i, "es par")) # Imprime i y que es par
} else { # Si no
print(paste(i, "es impar")) # Imprime i y que es impar
}
}
# Ejemplo con next:
for (i in 1:10){
if (i %% 2 == 1){ # Si el resto de dividir i entre 2 es 1
next # Siguiente
}
print(i) # Imprime i
}2.2.3 El uso de condicionales y bucles sobre bases de datos
Para realizar esta práctica vamos a utilizar los datos que usé para una investigación pasada que está publicada en este enlace. Preparamos el chunk de R para cargar los datos (recuerda que para cargar los datos en R hay que seleccionar una ruta determinada):
# Cargamos los datos
datos <- openxlsx::read.xlsx("data/Dataset.xlsx", sheet=2) # Otra forma de "invocar" funciones de paquetes sin pasar por library. Recuerda **adaptar tu ruta**
head(datos) # Para ver las primeras filas. Se puede configurar como head(x, 3), etc. ID age gender hiseilevel location cook degree dairy olivesnuts herbsspices
1 1 20 1 1 1 0 0 3.421918 0.7123288 0.27945205
2 2 45 1 2 1 1 0 6.643836 2.0000000 2.00000000
3 3 20 0 3 1 0 1 3.421918 0.4273973 0.00000000
4 4 18 1 3 1 0 0 1.498630 0.1315068 0.06575342
5 5 18 0 1 2 1 0 1.263014 0.0000000 1.21369863
6 6 18 1 2 1 0 0 9.627397 0.2794521 0.27945205
fruit veget oliveoil breadpasta potato redmeat sweets
1 4.115068 1.397260 0.49863014 2.0684932 1.495890 5.753425 6.3671233
2 4.279452 7.852055 2.50136986 5.1424658 0.000000 11.832877 1.9561644
3 1.463014 1.827397 2.50136986 5.2082192 0.460274 8.860274 18.8712329
4 1.895890 2.608219 0.49863014 1.6904110 0.460274 12.312329 8.8986301
5 3.610959 2.936986 0.06575342 0.9863014 0.460274 5.753425 0.9205479
6 5.117808 2.043836 1.00000000 2.0684932 1.495890 20.328767 4.2575342
whitemeat fish egg legumes alcoholicUA fastfood precooked
1 1.4958904 5.408219 1.495890 3.4520548 0.210 0.3287671 1.7643836
2 4.9863014 3.452055 4.986301 5.4082192 0.000 0.2630137 2.5671233
3 1.9561644 3.912329 1.495890 2.8767123 0.000 0.9205479 1.2547945
4 1.9561644 4.717808 3.490411 2.3013699 0.000 0.4602740 1.5260274
5 0.9205479 9.819178 0.460274 0.9205479 3.382 0.7068493 1.0246575
6 6.4821918 6.789041 3.490411 2.4164384 0.000 0.7616438 0.9041096
str(datos)'data.frame': 593 obs. of 24 variables:
$ ID : num 1 2 3 4 5 6 7 8 9 10 ...
$ age : num 20 45 20 18 18 18 18 18 18 18 ...
$ gender : num 1 1 0 1 0 1 0 0 1 0 ...
$ hiseilevel : num 1 2 3 3 1 2 3 1 2 2 ...
$ location : num 1 1 1 1 2 1 2 3 1 1 ...
$ cook : num 0 1 0 0 1 0 0 1 0 0 ...
$ degree : num 0 0 1 0 0 0 0 0 0 0 ...
$ dairy : num 3.42 6.64 3.42 1.5 1.26 ...
$ olivesnuts : num 0.712 2 0.427 0.132 0 ...
$ herbsspices: num 0.2795 2 0 0.0658 1.2137 ...
$ fruit : num 4.12 4.28 1.46 1.9 3.61 ...
$ veget : num 1.4 7.85 1.83 2.61 2.94 ...
$ oliveoil : num 0.4986 2.5014 2.5014 0.4986 0.0658 ...
$ breadpasta : num 2.068 5.142 5.208 1.69 0.986 ...
$ potato : num 1.5 0 0.46 0.46 0.46 ...
$ redmeat : num 5.75 11.83 8.86 12.31 5.75 ...
$ sweets : num 6.367 1.956 18.871 8.899 0.921 ...
$ whitemeat : num 1.496 4.986 1.956 1.956 0.921 ...
$ fish : num 5.41 3.45 3.91 4.72 9.82 ...
$ egg : num 1.5 4.99 1.5 3.49 0.46 ...
$ legumes : num 3.452 5.408 2.877 2.301 0.921 ...
$ alcoholicUA: num 0.21 0 0 0 3.38 ...
$ fastfood : num 0.329 0.263 0.921 0.46 0.707 ...
$ precooked : num 1.76 2.57 1.25 1.53 1.02 ...
Quizá se puede hacer de otra forma:
O de otra forma (de las más utilizadas actualmente):
datos <- read.xlsx("data/Dataset.xlsx", sheet=2)
library(dplyr)
datos2 <- datos %>% # Nuestro primer pipe :) -> Pongo datos2 por tener otra base distinta
mutate(across(c(gender, hiseilevel, location, cook, degree), as.factor))
str(datos2)'data.frame': 593 obs. of 24 variables:
$ ID : num 1 2 3 4 5 6 7 8 9 10 ...
$ age : num 20 45 20 18 18 18 18 18 18 18 ...
$ gender : Factor w/ 2 levels "0","1": 2 2 1 2 1 2 1 1 2 1 ...
$ hiseilevel : Factor w/ 3 levels "1","2","3": 1 2 3 3 1 2 3 1 2 2 ...
$ location : Factor w/ 3 levels "1","2","3": 1 1 1 1 2 1 2 3 1 1 ...
$ cook : Factor w/ 2 levels "0","1": 1 2 1 1 2 1 1 2 1 1 ...
$ degree : Factor w/ 2 levels "0","1": 1 1 2 1 1 1 1 1 1 1 ...
$ dairy : num 3.42 6.64 3.42 1.5 1.26 ...
$ olivesnuts : num 0.712 2 0.427 0.132 0 ...
$ herbsspices: num 0.2795 2 0 0.0658 1.2137 ...
$ fruit : num 4.12 4.28 1.46 1.9 3.61 ...
$ veget : num 1.4 7.85 1.83 2.61 2.94 ...
$ oliveoil : num 0.4986 2.5014 2.5014 0.4986 0.0658 ...
$ breadpasta : num 2.068 5.142 5.208 1.69 0.986 ...
$ potato : num 1.5 0 0.46 0.46 0.46 ...
$ redmeat : num 5.75 11.83 8.86 12.31 5.75 ...
$ sweets : num 6.367 1.956 18.871 8.899 0.921 ...
$ whitemeat : num 1.496 4.986 1.956 1.956 0.921 ...
$ fish : num 5.41 3.45 3.91 4.72 9.82 ...
$ egg : num 1.5 4.99 1.5 3.49 0.46 ...
$ legumes : num 3.452 5.408 2.877 2.301 0.921 ...
$ alcoholicUA: num 0.21 0 0 0 3.38 ...
$ fastfood : num 0.329 0.263 0.921 0.46 0.707 ...
$ precooked : num 1.76 2.57 1.25 1.53 1.02 ...
Ahora, si queremos categorizar variables continuas en función de distintos niveles o cortes, podemos hacer lo siguiente:
# Utilizando condicionales if/else como en Excel:
datos <- openxlsx::read.xlsx("data/Dataset.xlsx", sheet=2)
datos$age_cat <- NA # Se inicializa la columna
for (i in 1:nrow(datos)) {
if (datos$age[i] < 20) {
datos$age_cat[i] <- "Muy joven"
} else if (datos$age[i] >= 20 & datos$age[i] < 25) {
datos$age_cat[i] <- "Menos joven"
} else {
datos$age_cat[i] <- "Ojo, que si haces deporte te lesionas"
}
}
# Utilizando una función un poco más potente
datos <- openxlsx::read.xlsx("data/Dataset.xlsx", sheet=2)
datos$age_cat <- ifelse(datos$age < 20,
"Muy joven",
ifelse(datos$age >= 20 & datos$age < 25,
"Menos joven",
"Ojo, que si haces deporte te lesionas"))
# Utilizando funciones de dplyr
library(dplyr)
datos <- openxlsx::read.xlsx("data/Dataset.xlsx", sheet=2)
datos <- datos %>%
mutate(age_cat = (if_else(
age < 20,
"Muy joven",
if_else(age < 25, "Menos joven", "Ojo, que si haces deporte te lesionas")
)))# Esto arroja un problema: solución
summary(datos$age_cat) # Es una variable de caracteres Length Class Mode
593 character character
library(dplyr)
datos <- openxlsx::read.xlsx("data/Dataset.xlsx", sheet=2)
datos <- datos %>%
mutate(age_cat = factor(
if_else(
age < 20,
"Muy joven",
if_else(
age < 25,
"Menos joven",
"Ojo, que si haces deporte te lesionas")
)
))
# Otra forma
datos <- openxlsx::read.xlsx("data/Dataset.xlsx", sheet=2)
datos <- datos %>%
mutate(age_cat = factor(case_when(
age < 20 ~ "Muy joven",
age >= 20 & age < 25 ~ "Menos joven",
age >= 25 ~ "Ojo, te he dicho"
)))3. Tratamiento de bases de datos, análisis exploratorio y visualización: del tratamiento básico al {tidyverse} y a {ggplot2}
3.1 Tratamiento básico de bases de datos
Sin quererlo, antes ya hemos hecho algún tratamiento de datos. Incluso ya hemos utilizado parte del tidyverse. Ahora vamos a empezar a tratar los datos para conocerlos un poco mejor y usar algunas funciones básicas.
Ahora vamos a trabajar con un elemento (que creo) que es relevante y que nos permite trabajar con datos que estén en formato largo (long o narrow) o ancho (wide). En función del tipo de análisis que queramos hacer a veces tendremos que usar los datos en un formato o en otro. En el libro de The Epidemiologist R Handbook denominan a este proceso pivotar y recomiendan ver este capítulo de R for Data Science.
library(pacman)
pacman::p_load(rio, # Importar archivos
here, # Localizar archivos
kableExtra, # Crear y manipular tablas complejas
tidyverse) # Paquete que incluye ggplot2, dplyr, tidyr, etc. Gestión de datos, como sabemos
# Importamos datos: https://github.com/appliedepi/epirhandbook_eng/raw/master/data/malaria_facility_count_data.rds
library(DT)
count_data <- import("malaria_facility_count_data.rds")
datatable(count_data, options = list(pageLength = 5, filter = "top", lengthMenu = c(5, 10, 50, 100)))¿Qué observamos? ¿Cómo son estos datos?
Estos datos podrían ser entradas de información desde distintas instalaciones sanitarias (columna newid como si fuese el número del informe). Es decir, cada fila es un “caso” (sujeto, observación para un ID) que incluye toda la información de este. A este tipo de datos se le considera en formato wide o ancho.
Por ejemplo, los datos que hemos trabajado antes sobre factores sociales y consumo de alimentos están en formato ancho.
Esta base de datos es, cito literal: "Each observation in this dataset refers to the malaria counts at one of 65 facilities on a given date, ranging from count_data$data_date %>% min() to count_data$data_date %>% max(). These facilities are located in one Province (North) and four Districts (Spring, Bolo, Dingo, and Barnard). The dataset provides the overall counts of malaria, as well as age-specific counts in each of three age groups - <4 years, 5-14 years, and 15 years and older."
De ancho a largo (wide-to-long)
# Pasamos de ancho a largo considerando, por ejemplo, la edad y el total
df_long <- count_data %>%
pivot_longer(
cols = c(`malaria_rdt_0-4`, `malaria_rdt_5-14`, `malaria_rdt_15`, `malaria_tot`)
)
# Forma 1
df_long %>%
datatable(df_long, options = list(pageLength = 5, filter = "top", lengthMenu = c(5, 10, 50, 100)))# Forma 2 con tydyselect helper
count_data %>%
pivot_longer(
cols = starts_with("malaria_")
)# A tibble: 12,152 × 8
location_name data_date submitted_date Province District newid name value
<chr> <date> <date> <chr> <chr> <int> <chr> <int>
1 Facility 1 2020-08-11 2020-08-12 North Spring 1 malari… 11
2 Facility 1 2020-08-11 2020-08-12 North Spring 1 malari… 12
3 Facility 1 2020-08-11 2020-08-12 North Spring 1 malari… 23
4 Facility 1 2020-08-11 2020-08-12 North Spring 1 malari… 46
5 Facility 2 2020-08-11 2020-08-12 North Bolo 2 malari… 11
6 Facility 2 2020-08-11 2020-08-12 North Bolo 2 malari… 10
7 Facility 2 2020-08-11 2020-08-12 North Bolo 2 malari… 5
8 Facility 2 2020-08-11 2020-08-12 North Bolo 2 malari… 26
9 Facility 3 2020-08-11 2020-08-12 North Dingo 3 malari… 8
10 Facility 3 2020-08-11 2020-08-12 North Dingo 3 malari… 5
# ℹ 12,142 more rows
# Forma 3 con posición de columnas (y hay otras formas)
count_data %>%
pivot_longer(
cols = 6:9
)# A tibble: 12,152 × 8
location_name data_date submitted_date Province District newid name value
<chr> <date> <date> <chr> <chr> <int> <chr> <int>
1 Facility 1 2020-08-11 2020-08-12 North Spring 1 malari… 11
2 Facility 1 2020-08-11 2020-08-12 North Spring 1 malari… 12
3 Facility 1 2020-08-11 2020-08-12 North Spring 1 malari… 23
4 Facility 1 2020-08-11 2020-08-12 North Spring 1 malari… 46
5 Facility 2 2020-08-11 2020-08-12 North Bolo 2 malari… 11
6 Facility 2 2020-08-11 2020-08-12 North Bolo 2 malari… 10
7 Facility 2 2020-08-11 2020-08-12 North Bolo 2 malari… 5
8 Facility 2 2020-08-11 2020-08-12 North Bolo 2 malari… 26
9 Facility 3 2020-08-11 2020-08-12 North Dingo 3 malari… 8
10 Facility 3 2020-08-11 2020-08-12 North Dingo 3 malari… 5
# ℹ 12,142 more rows
Nota que hemos pasado de 3.038 filas a 12.152 (ahora los datos son más largos). Concretamente, 4 veces más largos (3.038 *4 (3 grupo de edad + total), pero tienen menos columnas.
Este tipo de datos en formato largo se usan en ocasiones (depende de cada paquete) para hacer análisis de supervivencia. Los datos de panel también suelen estar en formato largo. En salud, podría ser el número de observaciones de un paciente en un hospital (ej: lo veo el lunes y veo cómo está, lo veo el martes y veo cómo está y cada fila representa un momento distinto)
De largo a ancho (long-to-wide)
Que el formato largo se relevante ocasiones, no indica que el formato ancho no lo sea. A veces, el formato ancho es más útil para hacer análisis. En el ejemplo de los factores sociales y demográficos y consumo de alimentos esa organización de datos nos permitió hacer modelos de regresión lineal múltiple.
Para transformar los datos de largo a ancho usaremos otros datos:
# Datos usados: https://github.com/appliedepi/epirhandbook_eng/raw/master/data/case_linelists/linelist_cleaned.rds
linelist <- import("linelist_cleaned.rds")
head(linelist) case_id generation date_infection date_onset date_hospitalisation
1 5fe599 4 2014-05-08 2014-05-13 2014-05-15
2 8689b7 4 <NA> 2014-05-13 2014-05-14
3 11f8ea 2 <NA> 2014-05-16 2014-05-18
4 b8812a 3 2014-05-04 2014-05-18 2014-05-20
5 893f25 3 2014-05-18 2014-05-21 2014-05-22
6 be99c8 3 2014-05-03 2014-05-22 2014-05-23
date_outcome outcome gender age age_unit age_years age_cat age_cat5
1 <NA> <NA> m 2 years 2 0-4 0-4
2 2014-05-18 Recover f 3 years 3 0-4 0-4
3 2014-05-30 Recover m 56 years 56 50-69 55-59
4 <NA> <NA> f 18 years 18 15-19 15-19
5 2014-05-29 Recover m 3 years 3 0-4 0-4
6 2014-05-24 Recover f 16 years 16 15-19 15-19
hospital lon lat infector source wt_kg
1 Other -13.21574 8.468973 f547d6 other 27
2 Missing -13.21523 8.451719 <NA> <NA> 25
3 St. Mark's Maternity Hospital (SMMH) -13.21291 8.464817 <NA> <NA> 91
4 Port Hospital -13.23637 8.475476 f90f5f other 41
5 Military Hospital -13.22286 8.460824 11f8ea other 36
6 Port Hospital -13.22263 8.461831 aec8ec other 56
ht_cm ct_blood fever chills cough aches vomit temp time_admission bmi
1 48 22 no no yes no yes 36.8 <NA> 117.18750
2 59 22 <NA> <NA> <NA> <NA> <NA> 36.9 09:36 71.81844
3 238 21 <NA> <NA> <NA> <NA> <NA> 36.9 16:48 16.06525
4 135 23 no no no no no 36.8 11:22 22.49657
5 71 23 no no yes no yes 36.9 12:60 71.41440
6 116 21 no no yes no yes 37.6 14:13 41.61712
days_onset_hosp
1 2
2 1
3 2
4 2
5 1
6 1
age_cat gender n
1 0-4 f 640
2 0-4 m 416
3 0-4 <NA> 39
4 5-9 f 641
5 5-9 m 412
6 5-9 <NA> 42
7 10-14 f 518
8 10-14 m 383
9 10-14 <NA> 40
10 15-19 f 359
11 15-19 m 364
12 15-19 <NA> 20
13 20-29 f 468
14 20-29 m 575
15 20-29 <NA> 30
16 30-49 f 179
17 30-49 m 557
18 30-49 <NA> 18
19 50-69 f 2
20 50-69 m 91
21 50-69 <NA> 2
22 70+ m 5
23 70+ <NA> 1
24 <NA> <NA> 86
df_wide2 <- na.omit(df_wide) # Aprovecho y muestro cómo se eliminan na's3.2 Análisis exploratorio de datos con {tidyverse} (among others) y {ggplot2}
Para este apartado se sigue la información de un capítulo del libro Fundamentos de ciencias de datos con R
Funcionamiento elemental de {ggplot2}
El paquete ggplot2nos permite crear gráficos que, honestamente, son una pasada. Una de las principales fuentes de información sobre el uso de este paquete lo podéis encontrar en el libro de Wickham y sus distintas ediciones (mira cuáles son las diferencias entre ellas, por ejemplo, la segunda incluía aspectos de Data Management).
De un curso reciente de Storytelling con R en la UCLM, me quedé con una imagen que presenta resumen clave de cómo funciona el paquete ggplot2 (que adaptaron de aquí, recurso también interesante):
Los gráficos se hacen de menos a más y se insertan distintas capas que ofrecen personalización. Hagamos un ejemplo sencillo también extraído de los materiales del mencionado curso para ver cómo funcionan las distintas capas:
library(ggplot2) # Recuerda que estos dos paquetes están en {tidyverse}
library(dplyr) #cargamos los datos 'starwars'
library(ggthemes)
starwars <- starwars # redudante, pero para que se vea que se carga
# Paso 1: lienzo
ggplot(data = starwars)# Paso 3: geometría
ggplot(data = starwars,
mapping=aes(x=height, y=mass))+
geom_point() # Ojo outlier! ¿Quién es? Busca en la base de datos# Paso 4: quitar outlier y representar de nuevo hasta aquí
starwars2 <- filter(starwars, name != "Jabba Desilijic Tiure")
ggplot(data = starwars2,
mapping = aes(x = height, y = mass)) +
geom_point()# Paso 5: distintas configuraciones
ggplot(data = starwars2,
mapping = aes(x = height, y = mass, label = name)) +
geom_point() +
geom_text()ggplot(data = starwars2,
mapping = aes(x = height, y = mass, color = mass)) +
geom_point(size = 2) +
labs(
title = "Star Wars: Altura y peso de los personajes",
x = "Altura (cm)",
y = "Peso (kg)",
color = "Peso (kg)"
)ggplot(data = starwars2,
mapping = aes(x = height, y = mass, color = gender)) +
geom_point(size = 2) +
labs(
title = "Star Wars: Altura y peso de los personajes",
subtitle = "por género",
x = "Altura (cm)",
y = "Peso (kg)",
color = "Gender",
caption = "Fuente: {dplyr}")+
scale_colour_viridis_d()+ # Cambiamos la paleta de colores
geom_abline(intercept=lm(mass~height, data=starwars2)$coefficients[1], slope=lm(mass~height, data=starwars2)$coefficients[2], color="darkblue")+ # Añadimos una línea de regresión
theme_economist()# Ojo: todo esto se puede hacer exactamente igual usando {dplyr}:
starwars2 %>%
ggplot(aes(x = height, y = mass, color=gender))+
geom_point(size=2)+
labs(
title = "Star Wars: Altura y peso de los personajes",
subtitle = "por género",
x = "Altura (cm)",
y = "Peso (kg)",
color = "Gender",
caption = "Fuente: {dplyr}")+
scale_colour_viridis_d()+ # Cambiamos la paleta de colores
geom_abline(intercept=lm(mass~height, data=starwars2)$coefficients[1], slope=lm(mass~height, data=starwars2)$coefficients[2], color="darkblue")+ # Añadimos una línea de regresión
theme_economist()En ese curso también nos mostraron un recurso web bastante útil: la galería de gráficos con R. Este recurso proporciona los MUCHÍSIMOS tipos de gráficos que se pueden crear en función de los datos que queremos analizar, así como sus definiciones y el código básico para obtenerlos (emulando una vignette o, mejor, una cheatsheet como esta). Pero sin duda, si no tenemos claro qué tipo de gráfico podemos obtener el recurso más relevante es el de from data to viz. Es más, el recurso anterior bebe de este último. Otro recurso similar puede encontrarse aquí.
Con todo, no podemos olvidar la siguiente máxima. Cuando queremos hacer un gráfico lo primero que tenemos que pensar es a qué audiciencia nos estamos dirigiendo ya que esto condiciona qué tipo de gráfico podemos usar: no es lo mismo una audiencia académica que la población general o que un decisor político; no es lo mismo dar una conferencia divulgativa que una científica. ¡Adáptate al contexto!
Let’s practice!
Para realizar el análisis exploratorio de datos elemental vamos a utilizar los datos sobre consumo de alimentos y factores sociales. Vamos a responder a estas tres cuestiones:
- ¿Cuál es la distribución del consumo de alimentos?
- ¿Cuál es el porcentaje hombres y mujeres que estudian o ciencias sociales o ciencias de la salud?
- ¿Cuál es la correlación entre el consumo de los distintos grupos de alimentos?
El exploratorio básico podría ser:
# Cargamos los datos
datos <- openxlsx::read.xlsx("data/Dataset.xlsx", sheet=2)
# Una función útil es la función summary() o explore_all()
summary(datos) ID age gender hiseilevel location
Min. : 1 Min. :17.00 Min. :0.0000 Min. :1.000 Min. :1.000
1st Qu.:149 1st Qu.:18.00 1st Qu.:0.0000 1st Qu.:1.000 1st Qu.:1.000
Median :297 Median :19.00 Median :0.0000 Median :2.000 Median :1.000
Mean :297 Mean :20.21 Mean :0.4199 Mean :1.755 Mean :1.583
3rd Qu.:445 3rd Qu.:21.00 3rd Qu.:1.0000 3rd Qu.:2.000 3rd Qu.:2.000
Max. :593 Max. :48.00 Max. :1.0000 Max. :3.000 Max. :3.000
cook degree dairy olivesnuts
Min. :0.0000 Min. :0.0000 Min. : 0.000 Min. :0.00000
1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.: 1.630 1st Qu.:0.06575
Median :0.0000 Median :0.0000 Median : 2.501 Median :0.21370
Mean :0.2917 Mean :0.2395 Mean : 2.933 Mean :0.34373
3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.: 3.707 3rd Qu.:0.42740
Max. :1.0000 Max. :1.0000 Max. :14.507 Max. :3.21370
herbsspices fruit veget oliveoil
Min. :0.00000 Min. : 0.000 Min. : 0.0000 Min. :0.0000
1st Qu.:0.06575 1st Qu.: 1.321 1st Qu.: 0.9534 1st Qu.:0.4986
Median :0.27945 Median : 2.395 Median : 1.7068 Median :1.0000
Mean :0.57098 Mean : 2.895 Mean : 2.2268 Mean :1.1760
3rd Qu.:0.71233 3rd Qu.: 3.819 3rd Qu.: 2.7781 3rd Qu.:1.0000
Max. :7.00000 Max. :18.419 Max. :26.9945 Max. :6.0000
breadpasta potato redmeat sweets
Min. :0.000 Min. : 0.0000 Min. : 0.000 Min. : 0.000
1st Qu.:1.282 1st Qu.: 0.4603 1st Qu.: 7.479 1st Qu.: 2.301
Median :1.918 Median : 1.4959 Median :10.912 Median : 4.871
Mean :2.280 Mean : 1.3234 Mean :12.908 Mean : 6.385
3rd Qu.:2.978 3rd Qu.: 1.4959 3rd Qu.:15.841 3rd Qu.: 7.978
Max. :9.430 Max. :17.5096 Max. :80.989 Max. :42.460
whitemeat fish egg legumes
Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000
1st Qu.: 1.496 1st Qu.: 2.992 1st Qu.: 1.496 1st Qu.: 1.841
Median : 2.992 Median : 4.833 Median : 1.496 Median : 2.992
Mean : 3.461 Mean : 5.610 Mean : 2.790 Mean : 3.500
3rd Qu.: 3.951 3rd Qu.: 7.000 3rd Qu.: 3.490 3rd Qu.: 4.488
Max. :33.005 Max. :35.019 Max. :31.510 Max. :22.956
alcoholicUA fastfood precooked
Min. : 0.0000 Min. :0.0000 Min. : 0.0000
1st Qu.: 0.0000 1st Qu.:0.3288 1st Qu.: 0.5260
Median : 0.2100 Median :0.5589 Median : 0.9041
Mean : 0.7046 Mean :0.6404 Mean : 1.0286
3rd Qu.: 0.8400 3rd Qu.:0.9041 3rd Qu.: 1.3315
Max. :19.1250 Max. :2.8466 Max. :12.2630
datos %>% explore_all()datos |> explore_all() # Lo mismo que arribaexplore_all(datos) # Lo mismo que arriba# Otras estadísticas
mean(datos$legumes) # media de servings de legumbres[1] 3.499822
sd(datos$fastfood) # desviación típica de servings de fastfood[1] 0.4263786
# Como un summary más desarrollado y usando dplyr()
datos %>% select(dairy, fastfood, legumes) %>% descr()Descriptive Statistics
datos
N: 593
dairy fastfood legumes
----------------- -------- ---------- ---------
Mean 2.93 0.64 3.50
Std.Dev 2.02 0.43 2.75
Min 0.00 0.00 0.00
Q1 1.63 0.33 1.84
Median 2.50 0.56 2.99
Q3 3.71 0.90 4.49
Max 14.51 2.85 22.96
MAD 1.49 0.41 2.22
IQR 2.08 0.58 2.65
CV 0.69 0.67 0.79
Skewness 1.84 1.23 2.71
SE.Skewness 0.10 0.10 0.10
Kurtosis 5.28 2.77 13.69
N.Valid 593.00 593.00 593.00
Pct.Valid 100.00 100.00 100.00
Para responder a la primera cuestión podemos usar histogramas u otro tipo de gráficos:
# 1. Distribución consumo alimentos: podríamos hacerlo de distintas formas e incluso mostrar varios histogramas a la vez en un mismo cuadro
ggplot(datos, aes(x = dairy)) +
geom_histogram(binwidth = 0.5, fill = "#69b3a2", color = "white", alpha = 0.8) + # Personalización de colores
labs(
title = "Histograma de Lácteos",
x = "Valores",
y = "Frecuencia"
) +
theme_minimal(base_size = 14) + # Tema minimalista
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 16), # Centrar el título
axis.title = element_text(face = "bold"), # Ejes con texto en negrita
panel.grid.major = element_line(color = "gray85"), # Líneas de cuadrícula
panel.grid.minor = element_blank() # Eliminar líneas menores
)# También podemos plantear gráficos de cajas
ggplot(datos, aes(y = dairy)) +
geom_boxplot(fill = "orange")# Para trabajar con más de una variable hay que hacer los datos en formato largo, ¡ea!
datos_largos <- datos %>%
select(dairy, veget, legumes) %>%
pivot_longer(
cols = c(dairy, veget, legumes)
)
ggplot(datos_largos, aes(x=name, y = value)) +
geom_boxplot(fill = "orange")# Un poco feo y podríamos hacerlo más moderno y lustroso
ggplot(datos_largos, aes(x = name, y = value)) +
geom_violin(color = "darkblue") +
geom_boxplot(alpha = 0.5,
fill = "darkblue") +
labs(title = "Consumo de alimentos diario",
x = "Alimentos",
y = "Número de raciones medias") +
theme_minimal()Para estudiar el porcentaje de hombres y mujeres en los distintos tipos de estudios podemos estudiarlo a través de tablas de frecuencias elementales o usar algún gráfico de mosaico. Como siempre, se puede mejorar todo lo que queramos. Un ejemplo de esto:
# 2. Porcentaje de hombres/mujeres que estudian sociales/salud
table(datos$gender, datos$degree) # para tener una tabla sencilla de frecuencias
0 1
0 237 107
1 214 35
prop.table(table(datos$gender,datos$degree),margin=1)
0 1
0 0.6889535 0.3110465
1 0.8594378 0.1405622
# Los gráficos de mosaico son mu bonicos: uno básico
datos_fac = data.frame(gender=as.factor(datos$gender), degree=as.factor(datos$degree))
ggplot(data = datos_fac) +
geom_mosaic(aes(x = product(degree, gender),
fill = degree)) # Uno más completo:
# Instalar y cargar los paquetes necesarios
if (!require(ggplot2)) install.packages("ggplot2")
if (!require(ggmosaic)) install.packages("ggmosaic")
# Crear un conjunto de datos de ejemplo
library(ggmosaic)
set.seed(123)
# Gráfico de mosaico más "apañao"
ggplot(data = datos_fac) +
geom_mosaic(aes(x = product(degree, gender), fill = degree), alpha = 0.9, color = "white") +
scale_fill_viridis_d(option="D") +
labs(
title = "Distribución del tipo de estudios por género",
subtitle = "Así más bonico",
x = "Género",
y = "Proporción",
fill = "Tipo de estudios"
) +
theme_minimal(base_size = 14) + # Tema limpio y profesional
theme(
plot.title = element_text(hjust = 0.5, size = 18, face = "bold"), # Título centrado
plot.subtitle = element_text(hjust = 0.5, size = 14, face = "italic"), # Subtítulo
axis.text.x = element_text(angle = 45, hjust = 1), # Rotar etiquetas del eje X
legend.position = "right", # Leyenda a la derecha
legend.title = element_text(face = "bold"), # Título de la leyenda en negrita
panel.grid.major = element_blank(), # Sin cuadrícula principal
panel.grid.minor = element_blank() # Sin cuadrícula secundaria
)- El estudio de la correlación de los distintos grupos de alimentos se puede abordar de varias formas:
# Correlación básica
cor(datos$fish, datos$redmeat, method="pearson")[1] 0.2305928
cor.test(datos$fish, datos$redmeat, method="pearson") # Con significación
Pearson's product-moment correlation
data: datos$fish and datos$redmeat
t = 5.7611, df = 591, p-value = 1.346e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.1529162 0.3054377
sample estimates:
cor
0.2305928
# Matriz de correlaciones
datos_num <- datos[8:24]
mcor <- cor(datos_num) # matriz de correlaciones
round(cor(datos[8:24]),2) # con dos decimales dairy olivesnuts herbsspices fruit veget oliveoil breadpasta potato
dairy 1.00 -0.01 -0.04 -0.05 0.00 -0.01 0.08 -0.02
olivesnuts -0.01 1.00 0.13 0.15 0.11 0.07 0.07 0.06
herbsspices -0.04 0.13 1.00 0.07 0.46 0.18 0.05 0.19
fruit -0.05 0.15 0.07 1.00 0.20 0.12 0.01 0.02
veget 0.00 0.11 0.46 0.20 1.00 0.17 0.00 0.13
oliveoil -0.01 0.07 0.18 0.12 0.17 1.00 0.13 0.05
breadpasta 0.08 0.07 0.05 0.01 0.00 0.13 1.00 0.04
potato -0.02 0.06 0.19 0.02 0.13 0.05 0.04 1.00
redmeat 0.12 0.12 0.08 -0.04 0.05 -0.01 0.09 0.02
sweets 0.01 0.06 0.03 0.00 -0.06 0.02 0.07 -0.01
whitemeat 0.06 -0.01 0.10 0.08 0.10 0.03 -0.01 0.08
fish 0.07 0.15 0.16 0.19 0.15 0.06 0.05 0.08
egg 0.07 0.07 0.20 0.03 0.15 0.11 0.00 0.15
legumes -0.05 0.14 0.12 0.21 0.18 -0.03 0.09 0.09
alcoholicUA 0.02 0.01 0.01 0.04 -0.05 -0.12 -0.06 -0.03
fastfood 0.00 0.04 -0.07 -0.12 -0.16 -0.11 0.09 -0.02
precooked 0.00 0.15 0.05 0.09 0.10 -0.01 0.02 -0.04
redmeat sweets whitemeat fish egg legumes alcoholicUA fastfood
dairy 0.12 0.01 0.06 0.07 0.07 -0.05 0.02 0.00
olivesnuts 0.12 0.06 -0.01 0.15 0.07 0.14 0.01 0.04
herbsspices 0.08 0.03 0.10 0.16 0.20 0.12 0.01 -0.07
fruit -0.04 0.00 0.08 0.19 0.03 0.21 0.04 -0.12
veget 0.05 -0.06 0.10 0.15 0.15 0.18 -0.05 -0.16
oliveoil -0.01 0.02 0.03 0.06 0.11 -0.03 -0.12 -0.11
breadpasta 0.09 0.07 -0.01 0.05 0.00 0.09 -0.06 0.09
potato 0.02 -0.01 0.08 0.08 0.15 0.09 -0.03 -0.02
redmeat 1.00 0.05 0.10 0.23 0.05 0.05 0.05 0.21
sweets 0.05 1.00 -0.12 -0.07 -0.03 -0.03 0.04 0.19
whitemeat 0.10 -0.12 1.00 0.23 0.21 0.08 0.04 -0.10
fish 0.23 -0.07 0.23 1.00 0.05 0.06 0.03 -0.05
egg 0.05 -0.03 0.21 0.05 1.00 0.05 0.03 -0.02
legumes 0.05 -0.03 0.08 0.06 0.05 1.00 -0.04 0.00
alcoholicUA 0.05 0.04 0.04 0.03 0.03 -0.04 1.00 0.09
fastfood 0.21 0.19 -0.10 -0.05 -0.02 0.00 0.09 1.00
precooked 0.05 0.08 0.05 0.12 0.02 0.05 0.01 0.17
precooked
dairy 0.00
olivesnuts 0.15
herbsspices 0.05
fruit 0.09
veget 0.10
oliveoil -0.01
breadpasta 0.02
potato -0.04
redmeat 0.05
sweets 0.08
whitemeat 0.05
fish 0.12
egg 0.02
legumes 0.05
alcoholicUA 0.01
fastfood 0.17
precooked 1.00
# Gráficos de correlación que pueden ser útiles para explorar los datos
corrplot(mcor, method="circle", diag=F) # se puede arreglar tanto como queramos# Este también es mu bonico, pero aquí tenemos muchos datos y habría que arreglarlo estéticamente
datos_num %>%
ggpairs(progress = FALSE)+
theme(plot.margin = unit(c(1, 1, 1, 1), "cm")) Para rematar un gráfico de barras que no tiene nada que ver:
4. Introducción a los análisis estadísticos aplicados a la Economía y al Análisis de Políticas Públicas
Para saber qué tipo de análisis debemos aplicar, es importante conocer los tipos de datos (tipo de variable) tenemos y cuáles se adaptan mejor a cada modelo. Con el tiempo y la experiencia, las transformaciones de datos son comunes. Esto no significa que tengamos que torturar los datos hasta que obtengamos el resultado que intuimos podría ser adecuado: al contrario. Antes de plantear cualquier tipo de modelo estadístico tenemos que tener muy claro el diseño de nuestra investigación y qué han hecho investigaciones con temáticas y/o diseños similares. Esto nos permitirá saber el estado de la cuestión, lo común y lo novedoso que podemos aportar. Las transformaciones de datos o la creación de variables proxies surgirán por distintas razones como, por ejemplo, alguna limitación en el diseño. (casi) Siempre que trabajemos con honestidad y transparencia, obtendremos resultados satisfactorios.
Con este cuadro resumen (preparado por una antigua y buena profesora que me dio estadística cuando cursaba el máster e imagino que inspirado en el libro de Bioestadística amigable, porque el libro lo conocí en el máster) podemos ver distintos tipos de análisis en función de nuestros datos:
Recuerda también que hay paquetes que te permiten importar a R directamente las bases de datos desde las webs y navegar desde el programa. Esto es tremendamente útil. Un ejemplo sería el paquete {eurostat}. Otro ejemplo, usado en finanzas, sería descargar datos desde Yahoo! Finance.
# Eurostat
library(eurostat)
eurostat::eu_countries# A tibble: 27 × 3
code name label
<chr> <chr> <chr>
1 BE Belgium Belgium
2 BG Bulgaria Bulgaria
3 CZ Czechia Czechia
4 DK Denmark Denmark
5 DE Germany Germany (until 1990 former territory of the FRG)
6 EE Estonia Estonia
7 IE Ireland Ireland
8 EL Greece Greece
9 ES Spain Spain
10 FR France France
# ℹ 17 more rows
# Yahoo! Finance
# Descarga de datos -------------------------------------------------------
inicio <- as.Date("2000-01-01") #fecha de inicio
final <- as.Date(Sys.Date()) #fecha de final
activos <- c("^IBEX") #activos que quiero cargar
getSymbols(activos,
source = "yahoo",
from = inicio,
to = final) #extrae los datos[1] "IBEX"
datos_ibex <- IBEX[, 2:4]
# Gráfico con ggplot2
library(ggplot2)
datos_ibex2 <- data.frame(fecha = index(datos_ibex), coredata(datos_ibex))
ggplot(datos_ibex2, aes(x = fecha)) +
geom_line(aes(y = IBEX.High, color = "Cierre"), size = 1) +
geom_line(aes(y = IBEX.Low, color = "Bajo"), size = 1) +
geom_line(aes(y = IBEX.Close, color = "Alto"), size = 1) +
labs(
title = "Evolución del IBEX35",
x = "Fecha",
y = "Precio",
color = "Leyenda"
) +
scale_color_manual(values = c(
"Cierre" = "darkblue",
"Bajo" = "green3",
"Alto" = "darkred"
)) +
theme_minimal() +
theme(
legend.position = "top",
legend.title = element_text(face = "bold"),
axis.text.x = element_text(angle = 45, hjust = 1)
)# Gráfico dinámico --------------------------------------------------------
# Esto toma las columnas de máx, mín y cierre
dygraph(datos_ibex,
xlab = "Periodo",
ylab = "Cotización",
main = "Cotización IBEX 35 del 2000 a la actualidad") %>% dyRangeSelector()Con todo, para trabajar algunos modelos estadísticos, usaremos varias bases de datos ya utilizadas.
4.0. Diferencias en medias, proporciones, etc.
Estos análisis son importantes y suelen incluirse en las tablas de descriptivos de cualquier trabajo académico. Nos dan información relevante y muchas veces ya nos hablan de su relevancia como factores asociados (estudios transversales) o relacionados (en estudios longitudinales) en modelos un poco más complejos.
t.test(datos$redmeat~datos$gender, var.equal=T) # t-Student. Si no son normales U Mann-Whitney
Two Sample t-test
data: datos$redmeat by datos$gender
t = -2.7066, df = 591, p-value = 0.006994
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-3.5164750 -0.5591071
sample estimates:
mean in group 0 mean in group 1
12.05248 14.09027
var.test(datos$redmeat~datos$gender) # Si son o no varianzas iguales
F test to compare two variances
data: datos$redmeat by datos$gender
F = 1.0018, num df = 343, denom df = 248, p-value = 0.9926
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.7930661 1.2600626
sample estimates:
ratio of variances
1.00184
chisq.test(datos$gender, datos$degree) # Para estudiar diferencias en frecuencias en tablas de contingencia
Pearson's Chi-squared test with Yates' continuity correction
data: datos$gender and datos$degree
X-squared = 22.126, df = 1, p-value = 2.554e-06
4.1 Regresión lineal y logística simple
Recuerda que antes de hacer cualquier regresión hay que estudiar cuáles son las condiciones necesarias para poder hacer incluir las variables en una regresión. En general, las regresiones se utilizan para predecir (mejor R^2), pero también es común usarlas para estudiar asociaciones (sin fines de predicción), donde el poder de explicación de la varianza no tiene por qué ser tan alto.
En la regresión lineal sus supuestos básicos son: - Linealidad - Independencia del error o residuos: no hay correlación entre errores - Normalidad de error: error normalmente distribuido - Homocedasticidad (ausencia de heterocedasticidad ;-) ): hace referencia a que la varianza de los errores de la regresión es la misma para cada observación (así la variabilidad de la respuesta será la misma para valores bajo y altos de la variable independiente)
# Lineal
model <- lm(datos$redmeat ~ datos$veget)
summary(model) # model es un objeto con info dentro! Abre con $
Call:
lm(formula = datos$redmeat ~ datos$veget)
Residuals:
Min 1Q Median 3Q Max
-14.552 -5.487 -1.793 2.823 67.569
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.5250 0.5100 24.557 <2e-16 ***
datos$veget 0.1721 0.1560 1.103 0.27
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.095 on 591 degrees of freedom
Multiple R-squared: 0.002055, Adjusted R-squared: 0.0003661
F-statistic: 1.217 on 1 and 591 DF, p-value: 0.2704
par(mfrow=c(2,3))
# Linealidad en residuals vs fitted, no hay patrones
# Independencia del error
library(lmtest)
dwtest(model) # Aprox a 2 no hay autocorrelación - Durbin Watson
Durbin-Watson test
data: model
DW = 1.9958, p-value = 0.4785
alternative hypothesis: true autocorrelation is greater than 0
# Normalidad del error: el gráfico Q-Q. PUntos cerca de la línea.
# Homocedasticidad: varianza de residuos constante. Con gráfico si resid vs fitted no tiene patrón - Si, por ejemplo forma de embudo, heterocedasticidad. El test: Breusch-Pagan o White
bptest(model) # BP - Sale heterocedasticidad - Se puede corregir aplicando transformaciones logarítmicas en las variables
studentized Breusch-Pagan test
data: model
BP = 4.361, df = 1, p-value = 0.03677
bptest(model, ~fitted(model)+I(fitted(model)^2)) # White: no asume relacion lineal entre residuos y variables independ.
studentized Breusch-Pagan test
data: model
BP = 4.3984, df = 2, p-value = 0.1109
Error in log(datos$redmeat ~ log(datos$veget)): Argumento no numérico para una función matemática
# Logística: esta la vemos por encimísima:
model2 <- glm(datos$degree~datos$sexo, family="binomial") # hay más familias, depende de los datosError in model.frame.default(formula = datos$degree ~ datos$sexo, drop.unused.levels = TRUE): tipo (NULL) inválido para la variable 'datos$sexo'
Da error, ¿por qué? Soluciónalo con lo que ya sabemos.
datos$degree_f <- as.factor(datos$degree)
datos$gender_f <- as.factor(datos$gender)
model2 <- glm(datos$degree_f~datos$gender_f, family="binomial")
model2 <- glm(degree_f~gender, data=datos, family="binomial")
summary(model2)
Call:
glm(formula = degree_f ~ gender, family = "binomial", data = datos)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.7952 0.1165 -6.828 8.62e-12 ***
gender -1.0154 0.2164 -4.693 2.69e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 652.84 on 592 degrees of freedom
Residual deviance: 628.70 on 591 degrees of freedom
AIC: 632.7
Number of Fisher Scoring iterations: 4
OR <- exp(coef(model2)) # OR y sus IC
IC <- exp(confint(model2))
data.frame(OR,IC) OR X2.5.. X97.5..
(Intercept) 0.4514768 0.3580582 0.5655336
gender 0.3622587 0.2344588 0.5486563
4.2 Regresión lineal múltiple
Para realizar la regresión lineal múltiple se utiliza información y bases de datos de un manual titulado Introductory Econometrics: A modern approach y las soluciones están aquí. El ejemplo sirve para preguntarse sobre el efecto de distintas variables en la cantidad de manzanas ecológicas demandadas por una familia, incluyéndose entre las variables independientes el precio de las manzanas ecológicas, el de las normales y otras de caracter socioeconómico. En este ejemplo se cumple la teoría básica de la demanda: efecto negativo del precio de la manzana eco y efecto positivo del precio de la manzana normal.
library(wooldridge)
apple <- wooldridge::apple
model_app <- lm(data=apple, ecolbs~ ecoprc+regprc+faminc+hhsize+educ+age)
summary(model_app)
Call:
lm(formula = ecolbs ~ ecoprc + regprc + faminc + hhsize + educ +
age, data = apple)
Residuals:
Min 1Q Median 3Q Max
-2.570 -1.146 -0.595 0.515 39.940
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.056770 0.892650 1.184 0.237
ecoprc -2.861237 0.591991 -4.833 1.68e-06 ***
regprc 3.006077 0.712308 4.220 2.79e-05 ***
faminc 0.002195 0.002865 0.766 0.444
hhsize 0.063094 0.067780 0.931 0.352
educ 0.034313 0.045314 0.757 0.449
age 0.001389 0.006763 0.205 0.837
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.486 on 653 degrees of freedom
Multiple R-squared: 0.04022, Adjusted R-squared: 0.0314
F-statistic: 4.561 on 6 and 653 DF, p-value: 0.0001519
# Sacado también de https://herioscarlanda.wordpress.com/wp-content/uploads/2018/10/wooldridge-2009-introduccic3b3n-a-la-econometrc3ada-un-enfoque-moderno.pdf4.3 Diferencias en diferencias (diff-in-diff)
El análisis diff-in-diff o DiD es útil para experimentos naturales o quasiexperimentales (randomizados es complicado). Este tipo de análisis nos permite estudiar diferencias entre un grupo de intervención y otro de control en función del tiempo, cuestiones geográficas o por valores tras la adopción de una política o tras la sucesión de algún suceso relevante. Permite comparación transversal (entre grupo de control e intervención) y antes-después (cómo estaba cada grupo antes-después del suceso). Es un método popular para hacer inferencia causal en microeconometría aplicada.
4.3.0. Introducción
Para estudiar este método usaré el desarrollo del modelo del manual Learning Microeconometrics with R de Adams (2021). Las medidas repetidas en un individuo permite infererir causalmente el efecto del tratamiento y descubrir características no observadas, siendo esto útil en conjuntos de datos de panel.
Considera que los datos son como los de la figura que aparece más abajo.
Estamos interesadxs en el efecto de X en Y que se denota por b. Este efecto varía por cada individuo (faltaría un sub “i”, vaya). Considera que observamos dos Xs diferentes para el mismo individuo: dos niveles de educación distinto. Así, podemos estimar el efecto de la educación del individuo en el resultado de interés (Y). Podemos medir la diferencia en la renta asociado a un aumento en la educación: el efecto causal, en otras palabras, de la educación sobre la renta. Sin embargo, si observamos el mismo individuo con dos niveles de educación distintos, estamos observando al individuo en dos veces diferentes. Esas caractéristicas no observadas (U) también pueden cambiar entre el tiempo: ha un efecto de la educación y del tiempo en Y. El efecto del tiempo se denota por d. Esta imagen muestra que la relación entre X e Y es confusa. Hay una backdoor relationship entre X e Y, que funciona a través del tiempo: puedes ir de X a Y o de X a tiempo, de tiempo a U y de U a Y. Si el efecto de Tiempo a X y a U es 1, la regresión de X a Y es b+d. Para estimar b necesitamos una forma de estimar d.
4.3.1 Primeras diferencias
Observamos un outcome de interés para el individuo i en dos momentos del tiempo, 1 y 2: \[ y_{it} = a + b x_{it} + d v_{it} + \epsilon_{it} \] donde el outcome está determinado por la variable política/intervención “x” y el efecto no observado del tiempo por “u”. Tomando primeras diferencias podemos considerar el efecto no observado del tiempo (u):
\[ y_{i2} - y_{i1} = b (x_{i2} - x_{i1}) + d (v_{i2} - v_{i1}) + \epsilon_{i2} - \epsilon_{i1} \]
Si asumimos que el efecto del tiempo es 1, se puede añadir un intercepto para medir el efecto del tiempo. Lo desarrollamos en el siguiente ejemplo
Ejemplo con datos de panel simulados y primera estimación
Considera estos datos de panel simulados:
Estos datos permiten observar el cambio de x para un individuo. Dado que estamos interesados en el efecto de este cambio en x podemos mirar primeras diferencias. Hacemos regresión del cambio en X sobre el cambio en Y.
Call:
lm(formula = diffy ~ diffx)
Residuals:
Min 1Q Median 3Q Max
-5.1334 -0.9315 0.0577 0.8596 4.7124
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.77099 0.09013 -41.84 <2e-16 ***
diffx 2.74810 0.15416 17.83 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.415 on 998 degrees of freedom
Multiple R-squared: 0.2415, Adjusted R-squared: 0.2408
F-statistic: 317.8 on 1 and 998 DF, p-value: < 2.2e-16
Si se usa el intercepto, b se aproxima a 3 y el intercepto (d) a -4.
4.3.2 Diff-in-diff
Ahora vamos a desarrollar un ejemplo para aplicar un diff-in-diff siguiendo los pasos proporcionados en este tutorial (el ejemplo también está desarrollado en el libro que estamos siguiendo).
Los datos del ejemplo de Card and Krueger (1994) consisten en unas encuestas en restaurantes antes y después de que el estado de New Jersey aumentase el salario mínimo. El objetivo era estudiar el impacto en el empleo. En 1991, era de $4.25/h y en 1992, de $5.05/h. Así, se podría comparar simplemente el impacto en el empleo de la variación del salario mínimo tras la aplicación normativa (antes-después). Sin embargo, no parece suficiente el antes-después sin más porque otros factores podrían haber intervenido. Por tanto, se utiliza como referencia el estado vecino de Pensilvania como grupo de control (o sin tratamiento), considerando que que los restaurantes de ese estado no sufrieron impacto por el cambio normativo en New Jersey. La teoría económica ortodoxa indica que subidas de salario mínimo conllevan pérdidas de empleo, ¿será así?
Empezamos explorando los datos ¿sabemos ya cómo hacerlo?:
viridis_color <- viridis(5)[1] # saco un color
summary(rest) bonus chain co_owned date
Length:820 Length:820 Length:820 Min. :33553
Class :character Class :character Class :character 1st Qu.:33919
Mode :character Mode :character Mode :character Median :33921
Mean :33924
3rd Qu.:33925
Max. :33969
NA's :410
empft emppt firstinc hrsopen
Min. : 0.000 Min. : 0.00 Min. :0.0000 Min. : 7.00
1st Qu.: 2.000 1st Qu.:11.00 1st Qu.:0.1500 1st Qu.:12.00
Median : 6.000 Median :17.00 Median :0.2000 Median :15.00
Mean : 8.239 Mean :18.75 Mean :0.2214 Mean :14.45
3rd Qu.:12.000 3rd Qu.:25.00 3rd Qu.:0.2500 3rd Qu.:16.00
Max. :60.000 Max. :60.00 Max. :1.0000 Max. :24.00
NA's :18 NA's :14 NA's :123 NA's :11
inctime meals ncalls nmgrs
Min. : 2.00 Length:820 Min. :0.000 Min. : 0.000
1st Qu.:13.00 Class :character 1st Qu.:0.000 1st Qu.: 3.000
Median :19.00 Mode :character Median :1.000 Median : 3.000
Mean :20.08 Mean :1.447 Mean : 3.452
3rd Qu.:26.00 3rd Qu.:2.000 3rd Qu.: 4.000
Max. :52.00 Max. :9.000 Max. :10.000
NA's :97 NA's :249 NA's :12
nregs nregs11 observation open
Min. :1.000 Min. :0.000 Length:820 Min. : 0.000
1st Qu.:3.000 1st Qu.:2.000 Class :character 1st Qu.: 6.500
Median :3.000 Median :3.000 Mode :character Median : 7.000
Mean :3.601 Mean :2.689 Mean : 8.076
3rd Qu.:4.000 3rd Qu.:3.000 3rd Qu.:10.500
Max. :8.000 Max. :7.000 Max. :11.500
NA's :28 NA's :39 NA's :11
pctaff pentree pfry psoda
Min. : 0.00 Min. :0.410 Min. :0.6700 Min. :0.410
1st Qu.: 15.00 1st Qu.:0.940 1st Qu.:0.8500 1st Qu.:1.000
Median : 50.00 Median :1.030 Median :0.9400 Median :1.050
Mean : 48.87 Mean :1.338 Mean :0.9315 Mean :1.045
3rd Qu.: 80.00 3rd Qu.:1.843 3rd Qu.:1.0100 3rd Qu.:1.090
Max. :100.00 Max. :3.950 Max. :1.3700 Max. :1.490
NA's :454 NA's :36 NA's :45 NA's :30
region sheet special state
Length:820 Min. : 1.0 Length:820 Length:820
Class :character 1st Qu.:119.0 Class :character Class :character
Mode :character Median :237.5 Mode :character Mode :character
Mean :246.5
3rd Qu.:372.0
Max. :522.0
status type wage_st emptot
Length:820 Length:820 Min. :4.250 Min. : 0.00
Class :character Class :character 1st Qu.:4.500 1st Qu.:14.50
Mode :character Mode :character Median :5.000 Median :20.00
Mean :4.806 Mean :21.03
3rd Qu.:5.050 3rd Qu.:25.50
Max. :6.250 Max. :85.00
NA's :41 NA's :26
pct_fte
Min. : 0.00
1st Qu.:15.38
Median :32.65
Mean :34.03
3rd Qu.:54.05
Max. :90.36
NA's :32
rest %>%
explore_all()rest <- rest %>% filter(!is.na(rest$wage_st))
tapply(rest$wage_st, list(rest$state, rest$observation), mean) February 1992 November 1992
New Jersey 4.612134 5.080849
Pennsylvania 4.630132 4.617465
rest_before <- subset(rest, observation=="February 1992") # Saco la información de febrero de 1992 y muestro la función subset()
rest_after <- subset(rest, observation=="November 1992") # Saco la información de noviembre de 1992
rest_before %>% ggplot(aes(x=wage_st))+ # Histograma de salario
geom_histogram(fill=viridis_color, color="black", bins=20)+
labs(title="Salario en ambos estados Feb 1992", x="Salario", y="Frecuencia")+
xlim(3,7)+
#ylim(0,350)+
geom_vline(xintercept=c(4.25,5.05), color="black")rest_after %>% ggplot(aes(x=wage_st))+ # Histograma de salario
geom_histogram(fill=viridis_color, color="black", bins=20)+
labs(title="Salario en ambos estados Nov 1992", x="Salario", y="Frecuencia")+
xlim(3,7)+
# ylim(0,350)+
geom_vline(xintercept=c(4.25,5.05), color="black")Una vez explorados los datos, procedemos a estudiar las diferencias de medias ode empleo entre grupos (primeras diferencias:
differences <- rest %>%
group_by(observation, state) %>%
summarise(emptot = mean(emptot, na.rm = TRUE))
differences# A tibble: 4 × 3
# Groups: observation [2]
observation state emptot
<chr> <chr> <dbl>
1 February 1992 New Jersey 20.4
2 February 1992 Pennsylvania 23.4
3 November 1992 New Jersey 21.5
4 November 1992 Pennsylvania 21.6
# Treatment group (NJ) before treatment
njfeb <- differences[1,3]
# Control group (PA) before treatment
pafeb <- differences[2,3]
# Treatment group (NJ) after treatment
njnov <- differences[3,3]
# Control group (PA) after treatment
panov <- differences[4,3]Ahora sacamosel Efecto Medio del Tratamiento (Average Treatment Effect, ATT):
(njnov-njfeb)-(panov-pafeb) # Diferencia entre la diferencia de noviembre y febrero en cada estado emptot
1 2.876178
(njnov-panov)-(njfeb-pafeb) # Diferencia entre la diferencia de los estados en noviembre y febrero emptot
1 2.876178
Por último, podemos obtener el estimador DiD:
# Arreglamos los datos
rest <- mutate(rest, time=ifelse(observation=="November 1992", 1,0),
treated=ifelse(state=="New Jersey", 1,0))
did_model <- lm(emptot ~ time + treated + time:treated, data = rest)
summary(did_model) # Nota que los resultados no coinciden con el tutorial porque eliminamos NAs antes
Call:
lm(formula = emptot ~ time + treated + time:treated, data = rest)
Residuals:
Min 1Q Median 3Q Max
-15.929 -6.419 -0.972 4.528 64.581
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.429 1.079 21.714 <2e-16 ***
time -1.823 1.542 -1.183 0.2374
treated -3.010 1.203 -2.503 0.0125 *
time:treated 2.876 1.714 1.678 0.0938 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.282 on 756 degrees of freedom
(19 observations deleted due to missingness)
Multiple R-squared: 0.00892, Adjusted R-squared: 0.004987
F-statistic: 2.268 on 3 and 756 DF, p-value: 0.07932
El coeficiente de la itneracción tiempo-tratamiento es el estimador DiD. El tratamiento (aumento de salario mínimo) lleva en media a un incrempento de empleo en Nueva Jersey de 2.88 FTE (que es el número de trabajadores a tiempo completo (incluyendo gestores) más 0.5 veces el número de trabajadores a tiempo parcial)
4.4 Otros ejemplos de técnicas estadísticas aplicadas a economía y políticas públicas: análisis de supervivencia y transiciones entre estados
Estos ejemplos os los enseño desde mi propio PC pues son propios y los datos no son públicos.
5. Reproducibilidad con Quarto
Quarto es un sistema de publicación técnico y científico de código abierto. Permite crear documentos reproducibles e interactivos para distintos lenguajes: Python, R, Julia u Observable JS. En otras palabras, multilenguaje. Es la siguiente generación de R Markdown en Posit (RStudio) con otras características adicionales. Permite exportar en distintos formatos, pero con htmlya hay bastante juego.
Es el programa con el que he desarrollado esta hoja de ruta y os mostraré lo elemental para que podáis crear contenido similar. Lo (que creo que es más) importante:
- Configurar el YAML: basta con copiarlo de otro más desarrollado y modificarlo, probando qué pasa. Lo básico es título, fecha y autoría.
- Las almohadillas (#) son importantes para dividir e indexar el contenido.
- Los bloques de código se abren con tres acentos graves y el lenguaje de programación y se cierran con otros tres acentos graves (estos acentos son el opuesto al acento/tilde en castellano).
- Configurar bien la carpeta donde tenemos el proyecto pues se convertirá en el directorio de trabajo y podremos acceder fácilmente. Evitar rutas excesivamente largas porque si no puede fallar al leer las rutas de imágenes o bases de datos.
- Es otro lenguaje que funciona bastante bien con IA generativa. Haz las preguntas adecuadas y testea.
- Si te animas, únelo con github y comparte lo que desarrolles :).
Declaración de uso de IA generativa (no sé si hace falta ponerlo, pero creo que es lo correcto)
Para realizar este documento se han utilizado ChatGPT y Copilot (integración con GitHub) para autogenerar partes de código, para detección de errores, para corregirlos y para limpiar el código.